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, <^ ■ Abstract 

In inhomogeneous environments, the correct expression of the diffusive fiux is not always given 
by the Fick's law F = —D'Vn. The most general hydrodynamic equation modelling diffusion is 

r~| ' indeed the Fokker-Planck Equation (FPE). The microscopic dynamics of each specific system may 

affect the form of the FPE, either establishing connections between the diffusion and the convection 

c/3 ! term, as well as providing supplementary terms. In particular, the Fick's form for the Diffusion 

Q^l Equation may arise only in consequence of a specific kind of microscopic dynamics. It is also shown 

c/2 ' 

Q ■ how, in the presence of sharp inhomogeneities, even the hydrodynamic FPE limit may becomes 

^,^1 inaccurate and mask some features of the true solution, as computed from the Master Equation. 
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Introduction. The fluid modelling of the time- and space evolution of quantities within 
complex environments, whose dynamics may only be treated on statistical grounds, is made 
using the diffusion equation (DE) dtu = Dd'^n. This phenomenological equation arises from 
two more fundamental equations: the continuity equation for n: dtU = —d^T, and the Pick's 
law (or Fourier's law) [1] F = —Dd^n, where x and t are the spatial coordinate and the 
time, respectively; the diffusivity D is a constant dependent from the medium. A pedagog- 
ical overview of Fick's (Fourier's) law and diffusion equation may be found in [2|. 
The postulate of homogeneity may hold just as a first-order approximation, whereas most 
systems must ultimately allow for some degree of non-uniformity. Almost unavoidably, 
therefore, one is faced with the question: how DE has to be generalized to such systems. 
The exact answer to this question is of relevance for a plethora of problems in practically 
any branch of natural sciences: from physics, to chemistry, geology, biology, social sciences, 

Heuristically, the difficulty related to the generalization of DE may be understood as follows: 
an inhomogeneous environment should make D position-dependent: D -^ D{x). There are, 
however, several choices for F that differ when D = D{x), but that collapse to the same 
identical form when D is constant. Therefore, the problem may be restated as: what is the 
correct generalization of Fick's law (provided that one exists) in inhomogeneous environ- 
ments. 

This subject appears repeatedly addressed in literature; however, it is difficult to find the 
explicit exposition of a general solution. In Van Kampen's book |3], it is argued that one 
cannot decide a priori what the correct form for F is, which rather depends upon the prop- 
erties of the problem studied. Landsberg (jj] and references therein), points out that, to 
some extent, it is a matter of convention, provided that supplementary (convective) terms 
are added suitably. In other terms, the definition of a diffusive and a convective flux is not 
univocal, only the total flux is. The paper providing the clearest intuitive insight and at the 
same time detailed calculations about what goes on in such situations is probably Schnitzer's 
S. We „.ention *o the pape. QQ, featun,« eo^ute. expe.taents and p.ese„t.„g fu. 
ther bibliography about this subject. Papers [8|, |9|, HOj feature analytical and experimental 
work, demonstrating that the straightforward generalization of Fick's law F = —D{x)dxn{x) 
cannot hold in all systems. 

In order to quantitatively address the issue, it is necessary to deal with a reasonably 



accurate modelling of the dynamics at the microscopic level: transport equations, thus, will 

emerge at the level of large length scales. The tool we adopt is provided by the Master 

Equation (ME): 

dn{x,t) n{x,t) f , / / i n^i^'^t) 



+ [dx'p{x-x',x')^^^4lV- (1) 



dt t{x) 

ME (II]) yields a coarse grained probabilistic description of a microscopic system driven by 
a Markov process, and can be visualized as the continuity equation for the passive scalar 
quantity n(x,t) (which, properly speaking, is a probability density) subject to transitions 
("jumps") modifying its state from x' to x, with probability p(x — x', x'), and at a rate 1/t{x) 
(see chapter 1 of [11]). Equation ([T]) contains virtually all the solutions of the transport prob- 
lem, once the functions p and r are given. On the other hand, it is often unpractical to deal 
directly with it, particularly in higher-dimensional problems. Therefore, and particularly 
if a clear-cut separation of scales exists in the problem studied, it is customary to take its 
long- wavelength limit, which washes out details at the finest scales and turns the integral 



equation ([T]) into a 



amous differential equation: the Fokker-Planck Equation (FPE) (see. 



e.g., chapter 9 of |l2l|): 

Within the ME formulation, all the physics is built into the functions p and r. In the 
passage from ME to FPE, p and r are packed into the diffusive and convective terms, D, U. 
Therefore, the analytical expression of D, U, ultimately relies on the constraints that the 
problem to be solved places on p, r. Is it possible, basing upon general considerations on the 
microscopic dynamics, to identify equivalent classes of systems, that is, systems that lead to 
the same qualitative form of the FPE ? As we shall show later, the initial question advanced 
in this Introduction is related to this point: the Pick's form of the diffusion equation is a 
particular limiting case of the FPE, that arises when the microscopic dynamics fulfils a given 
simmetry. 

The purpose of this paper is to provide a discussion about this topic. Furthermore, we will 
address the broader issue of the validity of the scale separation at the basis of the FPE. We 
will show that, whenever, this hypothesis is not fulfilled, additional terms to the FPE need 
to be considered. 

From the Master Equation to the Fokker-Planck Equation. The simplest way to pass from 
ME to FPE is by expressing the integrand in Eq. ([1]) in terms of the small parameter 



A = x' — X, which is of order the mean jumping length Lp-. 

p{x-x',x') p{-A,x + A) . , .. .„s 

T[x') t{x + A) 

and expanding around x in powers of A (Kramer- Moyal expansion). However, this step 

is justified provided that p,T,n, are not strongly varying functions of x over distances of 

order Lp. If we assume that ra is a smooth function of x, we may concentrate on the other 

quantity: h = p/r. A branching into two cases is possible: (1) h is a smooth function, or 

(2) h is not. Although, condition (2) actually contains (1) as a particular case, it turns out 

convenient to consider them separately, since (1) is easier to deal with. 

Finally, we will consider also the case (3), when n itself is not a smooth function. 

Case (1): both n and h are smooth functions. We are allowed to make a Taylor expansion 

in powers of the function h x n. The result, truncated to second order, yields Eq. ([2]) with 

P(A,a;)^ ^_ 1 f ^p(A,x)^2 



U = / dA "^^ // A , D = - dA "^^ .\ A (4) 

J r(x) 2 J r(x) ^ ^ 

imiting the truncation to second order is ordinarily justified on the basis of Pawula theorem 



i3i,y. 



All the information relevant to our problem is packed into U, D. Two important cases are 
(^4) U = (dD/dx), or (B) U = 0. Case (^4) recovers Pick's law, while case (B) yields the 
solution 

dtn = dl{D{x)n{x)) (5) 

Both results may be verified by direct substitution into Eq. ([2]). It turns out that relation 
(A) arises straightforwardly from ME ([1]) by postulating the simmetry 

p{x' — x,x) p{x — x',x') p{A,x) p(— A,a;-|-A) ,^ 

t{x) ~ t{x') ~^ t{x) ~ t{x + A) ^ ^ 

which ensures the time reversal symmetry of the microscopic dynamics. Indeed, a first-order 
Taylor expansion of the second argument around x yields, after rearranging, 

.dp(—A,x) .. s ,.,.,. ,(i In r ,„. 

A ^^ ' ^ =p{A,x)~p{-A,x) + Ap{-A,x)—— (7) 

Using Eq. ([7]) into the integral s (H I) yields the sought result (A) (For a different derivation, 
see Prof. Peder's lecture notes |l5|). 

The solution (B) has some relevance, too, since it corresponds to the choice of a symmetrical 
kernel: p{A,x) = p{—A,x). Although apparently natural, the range of validity of this 



condition is actually rather narrow, as it cannot hold under smoothly varying conditions, 
where p{A,x) ^ p{—A,x); that is, the probability for a particle of jumping rightwards or 
leftwards cannot be the same. In order to better understand this point, let us consider 
a system where test particles collide against some scattering centres. Jumps are arcs of 
ballistic motion between two collisions. If the system is not homogeneous the density n^c of 
the scattering centres is not uniform. Let us suppose, say, dngc/dx < 0: a test particle at x 
has a larger probability of striking a scattering centre that is on its left {x — 5x) rather than 
on its right (x + 5x) , and therefore of being backscattered in the opposite direction. Hence, 
there is a larger probability of bouncing back rightwards than the converse. 
Having ruled out the case (B) for several inhomogeneous systems, one could wonder how 
general is condition (A). It turns out that (A) generically holds for a large class of 1-degree- 



of-freedom Hamiltonian systems 
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18| (see also [13] for another particular case). For 
more general systems, and especially in systems with more degrees of freedom, the above 
constraints between U and D cannot be guaranteed to hold any longer, and FPE may allow 
in principle for a wide variety of cases. Two such instances, recalled in [l6|], are: the self- 
consistent motion of charged particles in a set of Langmuir waves, and the 2-dimensional 
guiding-centre motion of a particle in a varying stochastic electrostatic field. 
Case (2): n is a smooth function hut h is not. Let us consider case (2), when p and/or r 
present sharp variations: we mean they vary on scales smaller than Lp. This is the situation 
when one needs modelling systems characterized by sudden transitions between regions with 
widely different physical properties. Therefore one is forced to study the case when we can 
still expand n in powers, but now must leave h unexpanded: after some calculations 

^"^ ^' ^bn) - ^ ((In) + Wn (8) 



dt dx'^ \ J dx 



2 

U = rhi + dxm2 

W = ^^ + d^mi + rho 

m, = [dx'ix-xy^^^^^^^^-6,0^ (J = 0,l,2) 

We have now an "extended" Fokker-Planck Equation, due to the presence of an additional 
Wn term. One may feel unconfortable about the presence of this term, since at first sight 



it appears to spoil the conservation of matter: dN/dt = d J n/dt = 0, which is built into 
Eq. ([1]) and Eq. ([2]). By integrating Eq. ([8]) over x, instead, one has dN/dt = J Wn{x)dx, 
which is not granted a priori to be zero. We will show, instead, that it is exactly the case: 
the role of W is that of transferring matter from one point to another, rather than that of 
a net sink or source. Heuristically, it may be guessed on the basis of the fact that Eq. (^ 
is just an intermediate passage in the chain of calculations leading from Eq. ([1]) to Eq. ([2]). 
Since both the starting and the end expressions do conserve matter, also the intermediate 
step must. For simplicity, from here on, we will consider the case of constant r = 1; hence, 
only variations in p will be dealt with. Let xq be a point around which p shows a sharp 
variation. For the sake of simplicity, let us consider a step-like variation around Xq = 0: 

, M \pl{x-x') x'<0 
p[x — X , X j = < (9j 

I Pr{x — x') x' > 

Since, far from xq, the system is almost homogeneous, we may suppose that the jumping 
lengths are symmetrical: Pl,b.{.x — x') = PL,R(y\x — x'\). We show now that, under these 
hypotheses, W reverses sign around Xq = 0: W{—x) = —W{x). For brevity, we proceed 
to demonstrate only that rfiQ is an odd function of x. Calculations for the other two terms 
carry on along the same lines: 



mo(x) = / dx'p{x — x' , x') — 1 = / dx'pL^x — x') + / dx'pR^x — x') — 1 

J Jx'<0 Jx'>0 

= / dx'piix — x') + / dx'pji{x — x') — I dx'pL{x — x') — / dx'pji{x — x') — 1 

J J Jx'>0 Jx'<0 



1 + 1 — / dx'piix — x') — / dx'pnlx — x) 

Jx'>0 Jx'<0 

— I / dx'pL{x — x)+ I dxpR{x — x) — l 



(10) 



If in the last line we make the change of variables x —>■ — x, x' —>■ — x', we get that the term 
within parentheses is just mo(— x). The second fundamental property of W is that W is 
different from zero only over a region around x = of width a few jumping lengths: It is 
straightforward to show that rho = drhi/dx = d'^rh2/dx^ = when |x| >> Lp . On the other 
hand, n is just a weakly varying function over distances of order the jumping length: this 
is a postulate implicit in the Taylor expansion done when going from Eq. ([T]) to Eq. (JHl). 
Therefore, combining the above results, J W{x)n{x) dx ~ n(0) J W{x) dx, and j W dx = 0. 
This concludes the demonstration. 




FIG. 1: W from Eq. ([T2|), with gl = l,aR = 2.5. 

In order to be more quantitative, let us consider the paradigmatic case of a gaussian diffusion 
with sharp variations of the jumping length: 



p{x — x', x') 
cr(x) 



1 



exp 



v/27ra2(x') 

(Tl , a; < 
cT/j , a; > 

It is straightforward to compute W for this system: 



[X — X 



l\2 



2a2(a;') 



W = X — 
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(11) 



(12) 



Its profile is given in figure ([T]). 
We have available an excellent experimental test bench of this result: paper [1(| presents 
a study of tracer diffusion between gelatine solutions with different viscosity, which means 
different effective diffusivity. The width of the interface between the two solution is very 
small, and can be considered as zero for our purposes. Hence, the whole system may be 
modelled as two regions with different jumping lengths, just done above. It is apparent at 
this point that applying Eq. ([2]) at this problem is of dubious validity, since p is discontinuous 
at x' = in its second argument. Nevertheless, just as an exercise, we may formally evaluate 
f/, D, for this choice of p: 



1 I (7? , X < 

f/ = 0, D = -^ 

^ ' cr^ , X > 



(13) 




FIG. 2: (Color online) Solution of Eq. §^ (solid line) and Eq. §^ (dashed line) at t = 32. Starting 
profile is flat: n{t = 0) = 1. The contoured region roughly envelops the area where experimental 
points lie (adapted from Fig. 5 of Ref. [iQl]). The two solutions are computed for slightly different 
boundary conditions: Eq. ([8]) has been numerically solved imposing dn/dx = at a; = ±40. Eq. 
([U may be integrated analytically. Its solution is discontinuous in x = 0, but the flux d{Dn)/dx 
is continuous p^]. For this function, dn/dx -^ only asymptotically as |x| — > oo. 



thereby recovering Eq. (I5l). W appearing in ([8]) has been computed in Eq. ( TT2l) . and D^ U 
are analytically computable alike, although we do not provide their complete expression here 
for saving space. The numerical values we choose for aL,R are ai = a/2.4, aR = V5.8, and 
solve for the time evolution of n{x,t) starting from a flat profile: n{x,t = 0) = 1. In Fig. 
([2]) we show the spatial profile at a later time for the solution of Eq. ([8]) as well as of Eq. 
dS]). These profiles should be compared against their counterpart from experiments (Fig. 5 
in ref. [1Q|). It becomes apparent that the smooth transition around x = 0, that exists in 
real data, is completely masked by the use of Eq. (j5]), while is correctly recovered by Eq. 



Case (3): n is not a smooth function and the full Master Equation is needed. This case 
acquires relevance when variations in n occur over scales smaller than the jumping length: 
Lp > Ln = n/ {dn/dx). In unbounded not driven systems this criterion may be satisfied 
only transiently, starting from highly localized profiles. Left to itself, density relaxes so 
as to fulfil Ln > Lp. However, L„ may be bounded by imposing absorbing boundaries. 
Hence, L„ < L, where L is the system' size. Since absorbing boundaries imply loss of 
density from the system, in order to maintain a steady state it is necessary to add a source, 
which is parameterized by its spatial extension, and hence a further typical length, Lg. The 
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FIG. 3: (Color online) Density profiles n due to the source (dashed curve) for three values of 
Lp = 0.0125 (black curve), 0.125 (green curve), 1 (orange curve). Here, the container has size 
L = 1, and the source width is L^ = 0.025. Reflecting conditions are imposed at x = and 
absorbing boundaries at x = 1. 

usual ordering is Lp < Lg < L, while now we investigate the reversal of this ordering. For 
simplicity, we will limit to consider the steady state d /dt = 0. Eq. ([T]) with a source 5* can 
be solved in Fourier space for g{x) = n{x)/T{x): 

S{kL,) 



9ik) 



(14) 



1 -p{kLp) 

In ( TT^ we have made it explicit that p depends upon scale length Lp, and S upon L^. 
For kLp,kLs << 1, one recovers the solution of the diffusion equation if p is a smooth 
symmetrical function: 1 —p{k) ^ k"^, k ^ 0. (If p lacks mirror symmetry, a convective term 
appears). Hence, over very large spatial scales one does not expect any novel feature to 
arise. However, when kLg ~ 1 but kLp » 1 the denominator of (IT4l) is almost one: p must 
have finite support, thus p ^ 0,k -^ oo. Hence 



g{k) ^ S{k) 



(15) 



Close to the source, and over distances of order of the source's width, the density profile 
matches exactly that of S. This is an effect totally unpredictable within the FPE formula- 
tion, that instead would smear the density throughout the whole system. In order to check 
numerically this prediction, we have solved Eq. ([1]) for a given source profile and different 
values of Lp. The results, fully confirming our analytical estimates, are shown in Fig. ([3]). 
Conclusions. In summary, inhomogeneity may have dramatic effect on the modelling of 



the spreading of some quantity into a medium. It affects the writing of the equation of mo- 
tion, turning it into a problem that has unambiguous solutions, but generally not valid ones 
for all classes of systems. Two fundamental criteria are identified. From the one hand, the 
microscopic dynamics imposes constraints between the diffusive and convective coefficients 
of the FPE. On the other hand, the form of the FPE itself depends on the existence of scale 
separation between typical lengths existing in the system: when this criterion is fulfilled and 
transport scales are clearly smaller than any other scale, the classical FPE is valid. On the 
opposite side, when no clear separation can be made, additional terms need to be added into 
an "augmented" FPE, till to the extreme case, where only the full ME may provide correct 
results. The first criterion is obviously the most important: FPE is often used for modelling 
purposes starting from a limited knowledge of the underlying microscopic dynamics. Clearly, 
one cannot establish a priori if scale separation does hold in the problem at hand. In these 
situations, therefore, the standard FPE has to be used, the only risk being that of missing 
some small-scale features. 
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